#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time    : 2018/1/31 23:25
# @Author  : Deyu Tian
# @Site    :
# @File    : config.py
# @Software: PyCharm Community Edition

import numpy as np
from math import *
import config
from image_tools import *


def _build_illum_model(slope, zenith, azimuth, aspect):
    """
    build the illumination terrain model
    :param slope:
    :param zenith:
    :param azimuth:
    :param aspect:
    :return:cos_i of incident angle
    """
    slope_rad = readSlopes(slope)
    demgt = readSlopecoeffs(slope)
    aspect_rad = readSlopes(aspect)
    print("max min of slope:", np.max(slope_rad), np.min(slope_rad))
    print("max min of aspect:", np.max(aspect_rad), np.min(aspect_rad))
    incident = np.zeros((slope_rad.shape[0], slope_rad.shape[1]), dtype=np.float32)
    #print("shape of slope:", slope_rad.shape)
    #print("shape of aspect:" ,aspect_rad.shape)
    for i in range(slope_rad.shape[0]):
        for j in range(slope_rad.shape[1]):
            if(slope_rad[i, j] < -5):
                #print("encountered nodata and killed...")
                incident[i, j] = -9999
                continue
            incident[i, j] = cos(slope_rad[i, j]) * cos(zenith) +\
                             sin(slope_rad[i, j]) * sin(zenith) * cos(azimuth-aspect_rad[i, j])
    print("max and min of cos_i:", np.max(incident), np.min(incident))
    np.save("{}\\incident.npy".format(config.basedir), incident)
    array2raster("{}\\incident.tif".format(config.basedir), demgt, incident)
    return incident




def cosine(surf, zenith, incident):
    """
    naive cosine method
    :param ToA:
    :param slope:
    :param aspect:
    :param zenith:
    :param azimuth:
    :return:
    """
    print("/////////////////////////begin running COSINE model///////////////////////////////////")
    cos_i = incident
    print("max and min of incident:", np.max(cos_i), np.min(cos_i))
    print("shape of incident", cos_i.shape)

    imggt, bands = readSixs(surf)

    print("bands number of ToA image:", len(bands))
    print("band shape:", bands[0].shape[0], bands[0].shape[1])
    for k in range(0, 6):
        surf_corr = np.zeros((cos_i.shape[0], cos_i.shape[1]))
        for i in range(cos_i.shape[0]):
            for j in range(cos_i.shape[1]):
                surf_corr[i, j] = bands[k][i, j] * cos(zenith) / cos_i[i, j]
        surf_corr[surf_corr > np.max(bands[k])] = np.max(bands[k])
        surf_corr[surf_corr < np.min(bands[k])] = np.min(bands[k])
        print("max and min of COSINE corrected 6S surf:", np.max(surf_corr), np.min(surf_corr))
        # np.save("{}\\ToA_corrected_B{}.npy".format(config.basedir, k+2), ToA_corr)
        array2raster("{}\\surf_corrected_B{}.tif".format(config.basedir, k + 2), imggt, surf_corr)
    print("/////////////////////end running COSINE model//////////////////////////////////////////")


def sun_canopy_sensor(surf, slope, zenith, indicent):
    """
    SCS method to generate modified ToA
    :param ToA:
    :param slope:
    :param zenith:
    :param incident:
    :return: modified ToA
    """
    cos_i = indicent
    print("/////////////////////begin running SCS model//////////////////////////////////////////")
    print("max and min of incident:", np.max(cos_i), np.min(cos_i))
    print("shape of incident", cos_i.shape)

    imggt, bands = readSixs(surf)
    slope_rad = readSlopes(slope)
    print("slope max and min:", np.max(slope_rad), np.min(slope_rad))

    print("bands number of ToA image:", len(bands))
    print("band shape:", bands[0].shape[0], bands[0].shape[1])
    for k in range(0, 6):
        surf_corr = np.zeros((slope_rad.shape[0], slope_rad.shape[1]))
        for i in range(slope_rad.shape[0]):
            for j in range(slope_rad.shape[1]):
                surf_corr[i, j] = bands[k][i, j] * (cos(slope_rad[i, j]) * cos(zenith)) / cos_i[i, j]
        surf_corr[surf_corr > np.max(bands[k])] = np.max(bands[k])
        surf_corr[surf_corr < np.min(bands[k])] = np.min(bands[k])
        print("max and min of SCS corrected 6S surf:", np.max(surf_corr), np.min(surf_corr))
        #np.save("{}\\ToA_corrected_B{}.npy".format(config.basedir, k+2), ToA_corr)
        array2raster("{}\\surf_corrected_B{}.tif".format(config.basedir, k+2), imggt, surf_corr)
    print("/////////////////////end running SCS model//////////////////////////////////////////")



def C_Factor(surf, zenith, indicent, c_values):
    """
    C-FACTOR method to generate modified surf
    :param ToA:
    :param slope:
    :param zenith:
    :param incident:
    :return: modified ToA
    """
    _C_values = c_values
    cos_i = indicent
    print("///////////////////////////begin running C-FACTOR model/////////////////////////////////////////////")

    print("max and min of incident:", np.max(cos_i), np.min(cos_i))
    print("shape of incident", cos_i.shape)

    imggt, bands = readSixs(surf)
    print("bands number of ToA image:", len(bands))
    print("band shape:", bands[0].shape[0], bands[0].shape[1])

    for k in range(0, 6):
        surf_corr = np.zeros((cos_i.shape[0], cos_i.shape[1]))
        print("max and min of band {}: {}, {}".format(k+2, np.max(bands[k]), np.min(bands[k])))
        for i in range(cos_i.shape[0]):
            for j in range(cos_i.shape[1]):
                surf_corr[i, j] = bands[k][i, j] * ((cos(zenith) + _C_values[k]) / (cos_i[i, j] + _C_values[k]))
        surf_corr[surf_corr > np.max(bands[k])] = np.max(bands[k])
        surf_corr[surf_corr < np.min(bands[k])] = np.min(bands[k])
        print("max and min of CC corrected 6S surf:", np.max(surf_corr), np.min(surf_corr))
        #np.save("{}\\ToA_corrected_B{}.npy".format(config.basedir, k+2), ToA_corr)
        array2raster("{}\\surf_corrected_B{}.tif".format(config.basedir, k+2), imggt, surf_corr)
    print("////////////////////////////////end running C-FACTOR model////////////////////////////////////////////////")



if __name__ == '__main__':
    zenith = 90 - 37.76268719
    azimuth = 150.08603472
    zenith_rad = zenith * pi / 180
    azimuth_rad = azimuth * pi / 180
    #cosine(config.surf6S, zenith_rad, config.indicent)
    #sun_canopy_sensor(config.surf6S, config.slope, zenith_rad, config.indicent)
    C_Factor(config.surf6S, zenith_rad, config.indicent, config._CC_regr_coeffs_and_intercepts)